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Abstract 

The derivation of the Feynman rules for lattice perturbation theory from ac- 
tions and operators is complicated, especially for highly improved actions such 
as HISQ. This task is, however, both important and particularly suitable for 
automation. We describe a suite of software to generate and evaluate Feyn- 
man rules for a wide range of lattice field theories with gluons and (relativistic 
and/or heavy) quarks. Our programs are capable of dealing with actions as 
complicated as (m)NRQCD and HISQ. Automated differentiation methods are 
used to calculate also the derivatives of Feynman diagrams. 
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LONG WRITE-UP 



1. Introduction 

Non-abelian gauge theories are the most important ingredient in our present 
understanding of elementary particles and their interactions. In particular, 
quantum chromodynamics (QCD) is now universally believed to be the cor- 
rect theory of the strong interactions. However, while perturbation theory has 
been used successfully in describing the scattering of particles by partons, the 
perturbative series does not converge at hadronic energy scales. Moreover, the 
phenomena of confinement and the hadronic spectrum arc fundamentally be- 
yond the reach of perturbation theory. Therefore, non-perturbative Monte Carlo 
simulations of lattice-regularised QCD are crucial in order to obtain a full de- 
scription and understanding of QCD phenomena. 

The lattice regularisation with a lattice spacing a does, however, introduce 
a sharp momentum cutoff at the momentum scale n/a. Connecting lattice mea- 
surements to their continuum counterparts therefore requires renormalisation 
factors accounting for the excluded high-frequency modes. In particular, renor- 
malisation is needed for QCD matrix elements, and for fixing the bare quark 
masses to be used in the lattice Lagrangian. Renormalisation is also necessary 
to determine the running of the strong coupling constant a s and to relate the 
lattice regularisation scale Ai at to Aqcd- Since the lattice regularisation also 
introduces discretisation errors, renormalisations arc also used to "improve" the 
lattice action in order to reduce these discretisation errors at a given lattice 
spacing. 

While in some cases the renormalisation constants can be determined non- 
perturbatively [TJI2], results at finite lattice spacing can depend upon the method 
used (cf. e.g. [3]), and non-perturbative methods do not cope well with oper- 
ators that mix under renormalisation. For these reasons, lattice perturbation 
theory plays an important role in determining the renormalisation constants 
needed to extract continuum predictions from lattice QCD. 

Given the breakdown of perturbation theory in low energy QCD, one might 
doubt whether it could work on the lattice. An argument in favour of its use 
is given in since the renormalisation factors may be thought of as compen- 
sating for the ultraviolet modes excluded by the lattice regulator, and since for 
typical lattice spacings a < 0.1 fm, the excluded modes have momenta in excess 
of 5 GeV, the running QCD coupling a s (7r/a) is small enough that perturbation 
theory should rapidly converge. The wide range of results reviewed in |S][S] show 
that perturbation theory is useful for a large range of lattice QCD processes. 
The assumption that non-perturbative effects do not contribute on these short 
length scales, can be tested directly in some cases by comparing higher order 
perturbative calculations with Monte Carlo simulations performed at a range 
of weak couplings [7J [HI M EH EH E2]), or by using so-called stochastic tech- 
niques [13 . In all these cases, the non-perturbative contributions turned out to 
be small. Other comparisons, such as [3], cannot distinguish non-perturbative 
effects from higher-order perturbative corrections. 
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Lattice perturbation theory therefore provides a reliable, and the only sys- 
tematically improvable, method for determining the full range of renormalisation 
constants [5]. 

As in the continuum, the calculation of lattice Feynman diagrams is a two- 
stage process. First, the lattice action and any operator insertions are Taylor- 
expanded in the (bare) coupling constant to give the propagators and vertices 
that form the Feynman rules (the "vertex expansion" stage). Secondly, these 
rules are then used to construct and numerically evaluate Feynman diagrams, 
possibly after some algebraic simplification (the "Feynman diagram evaluation" 
stage) . 

The latter task is more complicated than in the continuum due to the pres- 
ence of Lorcntz symmetry violating terms at finite lattice spacing, and on a 
finite lattice volume also by the more complicated nature of discrete momen- 
tum sums as compared to momentum integrals. Diagrams are therefore usually 
evaluated using numerical routines like Vegas [TH fTS], or proprietary math- 
ematical packages, possibly after manipulation using other computer algebra 
packages like Form [15] . 

The greater difficulty is posed, however, by the task of vertex expansion. 
Deriving the Feynman results on the lattice is far more complicated than in the 
continuum for a number of reasons. Firstly, lattice gauge fields are elements 
of the gauge Lie group itself rather than of its Lie algebra. To obtain the 
perturbative expansion of the action, we must therefore expand exponentials 
of non-commuting objects. As a consequence, the Feynman rules even for the 
simplest lattice action are already much more complicated than their continuum 
counterparts. 

Secondly, modern lattice theories generally contain a large number of addi- 
tional (renormalisation group irrelevant) terms chosen to improve specific as- 
pects of the Monte Carlo simulation, such as the rate of approach to the con- 
tinuum or chiral limits of QCD. Since there is no unique prescription for these 
terms, and the best choice depends on which quantities we are most interested 
in simulating, a large number of different lattice actions and operators are in 
use. Subtle though the differences between the lattice formulations may be, each 
choice provides a completely separate regularisation of QCD with its own set 
of renormalisation constants and in particular its own lattice Feynman rules. 
For a long time, the complications of the perturbative expansion have led to 
the calculations of renormalisation factors lagging far behind new developments 
in the improvement of lattice theories. In many cases this has restricted the 
physical predictions that could be obtained from the simulations. 

An automated method for deriving lattice Feynman rules for as wide a range 
of different theories as possible is therefore highly desirable. The vertex expan- 
sion should be fast enough not to impose undue constraints on the choice of 
action. To avoid human error, the user should be able to specify the action in a 
compact and intuitive manner. Since the evaluation of the Feynman diagrams 
can be computationally intensive, and will often need to be carried out in par- 
allel on costly supercomputing facilities, parsimony dictates that the Feynman 
rules should be calculated in advance on a different system, and rendered as ma- 
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chine readable files that can be copied to the supercomputer for the Feynman 
diagram evaluation stage. 

In this paper we describe a pair of software packages for deriving the Feyn- 
man rules for arbitrary lattice action^] and for evaluating the resulting vertices 
in a numerical Feynman diagram calculation. Our algorithm is based through 
our older algorithm [17 on the seminal work of Liischer and Weisz [TS]. A dif- 
ferent implementation of the latter has been used in [2U1 [H] , and a similar 
method is employed in [22 . 

The new feature of the algorithm presented here is that it is capable of 
expanding not only gluonic actions like the algorithm of 18J, and fermionic 
actions like our algorithm from |17j . but also far more complicated multiply- 
smeared fermionic actions with reunitarisation such as HISQ [23., and that it 
supports taking advantage of the factorisation inherent in some lattice actions, 
such as improved lattice formulations of NRQCD [24 . 

As in [IB] and [IT], the vertex expansion is performed completely indepen- 
dently of any boundary conditions, allowing for instance, the use of twisted 
periodic boundary conditions as a gauge-invariant infrared regulator |25] or 
for changing the discrete momentum spectrum in other ways [27J . 

We have used the software packages described for calculations of the renor- 
malised anisotropy in gauge theories [25J [S3] , to study the mean link in Landau 
gauge for tadpole improvement [TT], to measure the electromagnetic decays of 
heavy quark systems using NRQCD [3H1 [3H [3H 133] > to calculate the radiative 
corrections to the gluonic action due to Highly Improved Staggered (HISQ) sea 
quarks j34j [35] and the renormalisation of the self energy of heavy quarks using 
moving NRQCD (mNRQCD) [36]. 

The code is flexible and can be easily extended, as has already been done 
for lattice-regularised chiral perturbation theory [37] , perturbative calculations 
in the Schrodinger functional [31] EH] and anisotropy calculations [3D]. 

The structure of this paper is as follows. In the next Section, we review the 
basic expansion algorithm described in Ref. [17] . outlining some improvements 
in, for instance, the handling of automatic derivatives and spin matrices. In 
Sec. [3] we present novel extensions to the algorithm that are needed to describe 
complicated fermion actions, including HISQ, NRQCD and mNRQCD. 

Sec. [4]provides details of the implementation of the algorithm in the HiPPy 
and HPSRC codes and of their installation, testing and use. We make some 
concluding remarks in Sec. [5] Technical details are relegated to the appendices. 

1.1. Licence 

The HiPPy and HPSRC codes are released under the second version of the 
GNU General Public Licence (GPL v2). Therefore anyone is free to use or 
modify the code for their own calculations. 



x We shall use the term "action" so as to include measurement operators here and in the 
following. 
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As part of the licensing, we ask that any publications including results from 
the use of this code or of modifications of it cite Refs. [HI HZ] as well as this 
paper. 

Finally, we also ask that details of these publications, as well as of any bugs 
or required or useful improvements of this core code, would be communicated 
to us. 

2. Theoretical background 

In this section we describe the algorithms used to give the most efficient, 
yet generic, implementation of the Feynman rules for lattice actions. These 
extend the original work of Liischer and Weisz [T5] and developments described 
in Ref. [17]. 

2.1. Fields on the lattice 

We consider a D-dimcnsional hypercubic spacetime lattice with lattice spac- 
ing a and extent L^a in the /^-direction: 



where lattice sites arc labelled by a vector x € A. In the following, we will usu- 
ally set a = 1 (a lattice anisotropy can be introduced by rescaling the coupling 
constants in the action |28J). Let e M be a right-handed orthonormal basis, and 



A lattice path consisting of I links starting at site x can be specified by an 
ordered set of directions given by integers, S{ G {— D, . . . , — 1, 1, . . . , £>}: 



A 




,x D )eR D V/i6{l,. ..,£>}: ^e{0,...,^-l}} (1) 

a > 



y; a) = {x, y; s = [s , si, . . . , s/_i]} , 



(2) 



with the 7 



j point on the path being 




3 = 0, 
0<j<l, 



(3) 



and y = z\. 

For periodic boundary conditions, the momentum vectors are 




(4) 



and the Fourier expansion of a function <j> is 





X 



k 



where V = J\ is the lattice volume. 
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Twisted boundary conditions [H] provide a useful gauge-invariant infrared 
regulator in perturbative calculations [18*. These change the momentum spec- 
trum, converting colour factors into "twist matrices" associated with momenta 
interstitial to the reciprocal lattice (with the introduction of an additional quan- 
tum number, "smell" , for fermions [HI H21 [TTJ ) . The HPSRC code fully supports 
such boundary conditions but for simplicity we only discuss periodic boundary 
conditions in this paper. 

Following Ref. [18], we denote the gauge field associated with a link as 
Up >0 (x) G SU(N), and define U- M (x) = U}X x — ae^). The gauge potential 

e &\g(SU(N)) associated with the midpoint of the link is defined through 

W) = ex P (a 9 A» (x + £ c ,)) = £ ^ (6) 

r=0 

where g is the bare coupling constant. In terms of the anti-Hermitian generators 
of SU(N), 

A»=A«T a , [T a ,T b] = -f abc T c , K { T a T b ) = -\s ab . (7) 

Quark fields ip(x) transform according to the representation chosen for the gen- 
erators T a , which we take to be the fundamental representation (other choices 
will affect the colour factors, but not the structure of our algorithm). 

2.2. Perturbative expansion of Wilson lines 

The Wilson line L(x,y,U) associated with the lattice path C(x,y;s) is a 
product of links 

z-i ;-i 
L(x, y 7 U) ~Y\_ U H (zi) — II exp 

As all actions and operators can be written as sums of Wilson lines (possibly 
terminated by fermion fields) , our goal is to efficiently expand L in terms of the 
gauge potential in momentum space: 

L(x,y;A)=J2 ( ^f E "■ E ^(k,) . . . A« : (k r )x 

r fei iMi iOi k r ,fi r ,a r 

Vr(ki,fii,ai;...;k r ,fjbr,ar) • (9) 

The vertex functions V r factorise as 

V r {k 1 ,ni,a 1 ;...;k r ,fi r ,a r ) = C r {a\, . . . , a r ) Y r c (ki, fix; . . . ; k ri fj, r ) (10) 
with a momentum- and path-independent Clebsch-Gordan (colour) factor C r 

T 

C r ( ai ,...,a r )=l[T ai . (11) 

»=i 



n(si)agA\ Si \ (s 



(8) 
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It is therefore more efficient to calculate just the expansion of the reduced vertex 
functions, Y~ (with an appropriate description of the colour trace structure 
where ambiguous — see Appendix B of Ref. [17] for further details) . The reduced 
vertex function can be written as a sum of terms, each of which contains an 
exponential. For convenience, we will call each term a "monomial" : 

F r £ (fe 1 , / . 1 ;...;fe r , Mr )= E^^expfi^fe.-^f") , (12) 

n=l \ j=l / 

where for each combination of r Lorentz indices we have n r terms, each with an 
amplitude / and locations v of the r factors of the gauge potential, which are 
drawn from the locations of the midpoints of the links in the path C To avoid 
floating point ambiguities, we express the components of all position Z)-vectors 
as integer multiples of | (accounting for the factor of | in the exponent). 

In the HPSRC code, we use the convention that all momenta flow into the 
vertex, so Y%=i = 0- 

Eqn. ( |12[ ) makes clear that the number of monomials depends not just on 
the number of gluons r, but also on the choice of Lorentz indices {/^}, and that 
each monomial has a different amplitude and set of r positions. For clarity 
of presentation, we will, however, suppress these additional arguments in later 
expressions (notably Eqns. ([17] |T] |^j |44j [45]) ) . 

2.2.1. Implementation notes 

The generation of the Feynman rules for generic momenta thus reduces to a 
calculation of the amplitudes / and locations v of the monomials that build up 
the various reduced vertices Y r . 

This is all carried out in the HiPPy code, a description of which can be found 
in Sec. 4 of Ref. 1 1 7 | . The amplitudes and locations defining each monomial are 
encoded as an instance of the class Entity, and the collection of monomials that 
make up the reduced vertices is encoded as an instance of class Field. The data 
structures have been chosen to ensure that equivalent monomials are combined 
to minimise the size of the reduced vertex description. 

Once expanded, the monomials required for the reduced vertices at each 
order are written to disk as a text file. 

The HPSRC code reads these (previously generated) vertex files at runtime. 
For given momenta {k}, lorentz indices {/i} and colour indices, the Y r are 



constructed as given in Eqn. (12 1, which is then multiplied by the appropriate 



Clebsch-Gordan colour factor(s) to form the (Euclidean) Feynman rule, — V r . 

2.3. Realistic actions: the fermion sector 

Realistic lattice fermion and gauge actions require some refinements to this 
generic description. We begin with the fermion sector. The most general gauge- 
and translation-invariant action can be written as 

Spty, U)=J2J2 h ™ <P( X ) r w W(x, y, U) f(,(y) (13) 
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and consists of Wilson lines W denned by open paths W(x, y; s), each carry- 
ing an associated coupling constant /iyy and a spin matrix Tyy (possibly the 
identity) . 

Using the convention that all momenta flow into the vertex, the perturbative 
expansion is 

shim)=x;^ E ••• E A%(k 1 )...A%(k r )x 

E -0 fc (p) Vp,r(p, 9, c; fci,Aii,ai; fe r ,/i r ,a r )'0 c (q) . (14) 

p,q,b,c 

The Euclidean Feynman rule for the r-point gluon-fermion-anti-fermion vertex 
is —g r V Ft r, where the symmetrised vertex is: 

V Ftr (p,b;q,c;ki,Lii,ai; . . . ; k r , ju r , a r ) = 

■i E vC F , r {b,c;ai,...,a r ) cr ■Y F , r (p,q;k 1} fi 1 ; ...;k r ,fi r ) , (15) 

' a£S r 

where iS r is the permutation group of r objects and a € S r is applied to the 
gluonic variables, {fc}, {/i} and {a}. The normalisation factor of r! for this is 
additional to the r! factor arising from the Taylor expansion of the exponential 
in Eqn. ( 14 1 . The reduced vertex Yjr r = X/w ^w^fv is the sum of contributions 
from paths W. 

In most cases the Clebsch-Gordan colour factor is the matrix element: 

C F ,r{b,c;a x ,...,a r ) = (T 0l ...T 0r ) 6c , (16) 
and the reduced vertex function has the structure: 

Y F<r (p,q;kx,fJ,x;...;k r ,fi r )= E r n/n x 



exp (I ■ ^ + 9 ■ y + E fc j ' J J ' ( 17 ) 

where we understand T n = r r r^i n . Cases with more complicated colour struc- 
tures do arise, for example the use of traceless field strengths in QCD. Such 
structures are accommodated in the codes; monomials with different colour 
structures are distinguished using "pattern lists" (discussed in Appendix B of 
Ref . |17j ) and appropriate colour factors are applied to each when the Feynman 
rule is constructed. 

As there are no permutation symmetries in C F<r , there is no advantage to 
carrying out any symmetrisation in the HiPPy expansion code. In the HPSRC 
code, symmetrisation of the Feynman rule shown in Eqn. ( 15 1 carries a poten- 
tially significant computational overhead: the reduced vertices must be calcu- 
lated afresh for each permutation. Not all such permutations may be needed 
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because symmetries of a Feynman diagram can reduce the number of distinct 
contributions to its value from the terms in Eqn. (15). For this reason, sym- 
metrisation is not carried automatically in the HPSRC code and the user must 
therefore explicitly construct all permutations requiring a different calculation 
from Eqn. (15 1, applying the appropriate normalisation factor. 

2.4- Realistic actions: the gluon sector 
A typical gluonic action is 

x V 

built of Wilson loops P defined by closed paths V(x,x\ s), each with coupling 
constant c-p. The perturbative action is 



^) = Ett E •■• E A-(fci)...i£(*v)x 

t fei k r ,fi r .a r 

V Gtr (k ll nx, ai; . . . ; k r , /x r , a r ) . (19) 

The Euclidean Feynman rule for the r-point gluon vertex function is (— g r Vc!,r), 
and the vertex V G ,r is [IB] 

VG, r (kx,fj,x,ai; ...;k r ,fi r ,a r ) = 

"1 E CT ' C G,r{ai, ■■ - ,a r ) a - Y G , r {ki,Hi\ • ■ • ; Mr) , (20) 

' <reS r 

The reduced vertex = J^v c vYq r is the sum of contributions from paths V . 
As before, the (r!) factor normalises the symmetrisation. Yq t can be expanded 
as 

Y ciA k ^ Ml? • • • ; fe r> Mr) = E fn ex P I ^ E fcl ) ' ^ 21 ) 

n=l V i / 

In most cases we expect the lattice action to be real. Thus, for every monomial 
(fn, { v n,i}) m Eqn. plj ), there must be a corresponding term ((— l) r /*; {— w n ,t})- 
We can therefore speed up the evaluation of the Feynman rules by removing the 
latter term, and replacing the exponentiation in Eq. (21 1 with "cos" for r even, 
and with "i sin" for r odd. This can either be done by recognising conjugate 
contours in the action (e.g. S = \ Tr[P + P*]) and expanding only one, or by 
attaching a flag to each monomial to signal whether its complex conjugate has 
already been removed. 

If in addition the action has the form Eq. (18 1 with a single trace in the 
fundamental representation, the colour factors are 

C G ,r(*i> ■ ■ ■ > «»r) = \ I Tl ' ■ ■ ■ T a r ) + (~l) r Tr (T Qr . . . T Bl )] . (22) 
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which has a number of permutation symmetries: 



Cg,t = Xr{<?) C G ^ r , where XrO) 



1 for a a cyclic permutation, 

(— l) r for a the inversion. 

(23) 

There is thus a great advantage in carrying out some of the symmetrisation 



in Eqn. (20 1 at the expansion stage in the HiPPy code. Many of the extra 
monomials generated by symmetrising over subgroup Z r (generated by cyclic 
permutations and inversion) are equivalent and can be combined in the HiPPy 
code, significantly reducing the number of exponentiation operations required 
to construct the partially-symmetrised Y' G • 

VG,r(ki,fj,i,ai; ... ; k r ,{i r ,a r ) = ^ a ' c cA a ii —■> a r) x 

<j£S r /Z r 

a ■ F4 r (fei,/ii; k r ,fi r ) , (24) 

% = E ° ■ Ya, r ■ (25) 

•p 

o-£Z r 

The Xr(c) factors go into the amplitudes of the new monomials coming from 
the partial symmetrisation. 

The number of symmetrisation steps remaining to be carried out in the 
HPSRC code is the number of cosets in S r /Z r (one for r < 3, three for r = 4, 
twelve for r = 5 etc.). These symmetrisation steps (and the normalisation) are 
carried out automatically in the HPSRC gluon vertex modules. 

2.5. Diagram differentiation 

In many cases, such as when computing wavefunction renormalisation con- 
stants, one needs to calculate the derivative of a Feynman diagram with respect 
to one or more momenta. Whilst derivatives can be computed numerically 
using an appropriate local difference operator, such differencing schemes are 
frequently numerically unstable and require computing the Feynman diagram 
multiple times. Automatic differentiation methods [44] are a stable and cost- 
saving alternative. 



We can easily construct the differentiated Feynman vertex using Eqn. (12 1. 
If we want to differentiate with repsect to momentum component q v , we first 
construct a rank r object r = [t±, . .. ,r r ] which represents the proportion of 
momentum q in each leg of the Feynman diagram. Momentum conservation 
dictates Ti — 0. For instance, for a gluon 3-point function with incoming 
momenta (p, p + 2q, —2q), we would have r — [0, 2, —2]. The differentiated 
vertex is 



Y r c (ki,m; ...;fe n i«r)=^Y I ^TjV n j. tV J exp 



dq v 



n=l \i=l 




(26) 
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and so on for higher derivatives. We may therefore simultaneously calculate as 
many differentials as we need for the cost of just one exponentiation. If this 
momentum expansion is placed into an appropriate data structure for which 
appropriately overloaded operations have been defined, it is straightforward to 
create the Taylor series for a Feynman diagram by simply multiplying the vertex 
factors together. We use a slightly modified version of the TaylUR package 
[151 HB] to do this, which encodes the Taylor series expansion in the HPSRC 
Fortran code as a derived type taylor, for which all arithmetic operations and 
elementary functions have been overloaded so as to respect Leibniz's and Faa 
di Bruno's rules for higher derivatives of products and functions. 

In calculations that require only certain higher-order derivatives, the tay- 
lor multiplication can be significantly sped up by defining a mask that only 
propagates certain terms in the Taylor expansion. This has not, however, been 
implemented in the distributed version of the code. 

2.6. Spin algebra 

The Feynman rules for fermions contain spin matrices (which may be Pauli 
matrices, e.g. for NRQCD, or Dirac matrices for relativistic fermions). We can 
expand a generic spin matrix W using a basis {Ti}: 



The product of any two spin basis matrices Ti and Tj is another basis matrix 
(with label fiy) times a phase: 



We choose the basis so that these phases (f>ij are real. Where another convention 
is desired, the amplitudes Wi need to be adjusted appropriately. 

For Pauli matrices, we use a basis {Ti} — {l,iak} and I(W) C {0, . . . ,3}, 
giving: 



For the (Euclidean) Dirac matrices, we choose {Ti} = {1,7^,0^,/, 757^,75} and 
I(W) C {0, . . . , 15}. We define 75 = 71727374 and note the definition here 
Gp, v = \ [7^ , 7„] . The multiplication for the basic 7 matrices is thus 




(27) 



igJ(W) 



(28) 




(29) 




(30) 
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We can thus write the product of two general spin matrices as 



WW' 



E W 3 T 1 E W 'k T k 

K ]£I(W) j \k£l(W>) 

( \ 



E 

iGl(WW') 



E W j W 'k<l>jk 



(31) 



/ 



where I(WW) = {n l0 \ i £ I(W), j € I(W')}. 

We use this implicit representation of the spin matrices. Matrix opera- 
tions on explicit Dirac matrices take C(4 3 ) operations and introduce additional 
rounding errors. Eqn. (31 1, by contrast, needs |/(W)| x multiplications, 



depending on how many basis matrices are needed to describe W and W'. If 
|J(W)|, |J(W')| < 8, the latter method is more efficient and this is almost always 
the case. 

There are greater gains when inverting spin matrices of the form 



S 1 = a l + ^ °m7m j 



(32) 



as we might do when obtaining the propagator of a rclativistic quark (for 
NRQCD, the propagator is spin diagonal in the Pauli matrices and trivial to 
invert). The inverse is 



S = b l 



5> 



E«2 



(33) 



which is far more efficient than inverting a 4 x 4 matrix. Inversion of a gen- 



eral spin matrix (not of the form Eqn. ( 32 1 ) is less efficient with an implicit 
representation, but this is irrelevant in most perturbative calculations. 

Since all basis matrices except the identity are traceless, taking the trace of 
a spin matrix is a free operation in the implicit representation. 



2.6.1. Implementation notes 

In the HiPPy code, spin basis matrices are associated with monomials using 
an appropriate integer i, which is part of the Entity data structure. When terms 
are multiplied together, the factors <pjk are absorbed into the amplitude / of 
the resulting monomial. 

In the HPSRC code, we represent a spin matrix W as a defined type, spinor, 
which is encoded as a double array 

n 

(n;ii, . . .,i n ;wi, . . . ,w n ) = ^w k T ik . (34) 

fc=i 
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It turns out to be significantly more efficient for the order of terms in this array 
to not necessarily match a standard order of basis elements {1^}, hence the use 
of the index array ik- In particular, this allows us to omit basis elements with 
zero coefficient. 

Arithmetic operations have been overloaded to act appropriately on objects 
of this type, including implemention of the multiplication table. During inver- 
sion, an additional function argument, short_spinor, is used to employ the 
more efficient expression in Eqn. (33 1. 



3. Even more realistic fermionic actions 

Sec. [2] summarised the general method that was described in Ref. [T7]- In 
general, fermionic actions are much more complicated than gluonic actions and 
several algorithmic improvements are needed to efficiently calculate the associ- 
ated Feynman rules. We stress that by efficient we mean speed-ups of at least 
an order of magnitude. 

The algorithms described in this section can be used independently or to- 
gether, with the choice configurable at runtime of the code. All of these features 



have also been implemented using taylor- valued variables (as per Sec. 2.5 1 to 
provide automatic differentiation of Feynman rules. 



3.1. Summands and factors 

In many cases an action has a block- like structure that we can exploit to make 
the evaluation of the reduced vertices more efficient. This is particularly useful in 
the case of NRQCD and mNRQCD actions, which can be heuristically written 
as i^(l — ABCA)%j). Such an action can be expanded directly in the HiPPy 
code but the extremely large number of monomials makes this is inefficient (or 
impossibly slow and memory-hungry). It is clear why: there is often little scope 
for monomial compression between blocks. For instance, in (m)NRQCD, blocks 
AB and CA live on different timeslices of the lattice, and no compression is 
possible when combining them. 

Instead, we recognise that the blocks are combined in a gauge covariant 
manner, so that in AB, for instance, each contour in B starts where a contour 
in A finishes. Summing over the start/end location of each block we obtain 
a convolution of terms and can use the convolution theorem to construct the 
overall reduced vertex (i.e. the momentum-space Fourier transform) from those 
of the individual blocks. 

We refer to the action as being a sum over terms that we call "summands" . 
In the above example, there are two. Each summand is the convolution of 
a number of "factors", with one factor in the first summand and four in the 
second. 

The overall reduced vertex Yp r is the sum of the reduced vertices for each 
summand. For each summand, Yp. r is calculated by combining the reduced 

(k) 

vertices Y F ' for each of N factors that make up that summand, k = 1 . . . N. 
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Table 1: The elements P in the set of ordered partitions, P' 1 "', of the ordered set of the first 
r integers, {1,2, ... ,r}, for r = 1, 2 and 3. 



r P 


1^1 


r P 


1^1 


1 {{1}} 


1 


3 {{1,2,3}} 


1 


2 {{1,2}} 


1 


{{1,2}, {3}} 


2 


{{1},{2}} 


2 


{{1},{2,3}} 


2 






{{1},{2},{3}} 


3 



We generate these by expanding each factor of the action separately in the 
HiPPy code, with the convolution then carried out in the HPSRC code. 

Here we give the method for constructing Yp, r for a summand with N factors 
for general r. Expressions for specific r = . . . 3 are given in Appendix |A") 

In giving the general expression, we first establish some useful notation. 
Consider the ordered set of the first r integers: {1,2, ... ,r}. We can form an 
ordered partition of this set: {Px,P%, . . . , P z }, where the cardinality z = |{...}| is 
the number of elements in the partition. The set of all such partitions we denote 
p( r K For instance, the partitions P( r ' for r = 1, 2, 3 are shown in Table [I] Note 
that we do not consider unordered partitions (e.g. {{1,3}, {2}}) because the 
gauge fields are explicitly ordered in the paths (Wilson lines) making up the 
action. 

The general reduced vertex for a summand with N factors is then: 



Y F:r (p,q;ki,ni; . . .;k r ,ii r ) = 

e e (n 

PeP(r) l<n 1 <n 2 <---<n ]P] <N [QeP 

X Y F\\Q\ (P« ' -Pi+l i k<9i , MQi i 



n *#oWpo 



; kQ\Q\ ) ^Q\Q\ ) 
N 

n ^(-9.9) 



(35) 



where i = 1 . . . \P\ is the position of element Q in the ordered set P and no = 0. 
Momentum is conserved in the diagram, so p\ — p and subsequent Pi+\ = 
Pi + k,Q (with i defined as above) . Here k,Q refers to the summed momenta for 
the set Q (e.g. fe{i,2} = fei + ^2), implying = -q. 



5.^. Two-level actions 

Fermion actions often use fattened links to reduce discretisation errors in 
numerical simulations. By fattening or smearing a link, we will, in general, end 
up with an N x N complex colour matrix M, which can be expanded in powers 
of the gauge field A M . It is often the case that this matrix is reunitarised by 
projecting it back onto the gauge group SU(N) or, more usually, simply back 
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onto the related group U(N). Fattened links can then be further fattened, in 
an iterative procedure. An example of this is the HISQ action. 

To complement these numerical simulations, we need to do perturbative 
calculations using the same actions. We confine our attention here to the HISQ 
action (and simpler variants of the same form, for testing), with an iterated, 
two-level smearing procedure with an intermediate reunitarisation: 

U msQ = (F ASQ , o p m o F FAT7 ) [U] (36) 

where U = exp(gA) is the unsmeared gauge field, Pu(3) denotes the polar pro- 
jection onto U(3) (as used in simulations, and not 5*17(3) |47]), and the FAT7 
and ASQ' (a slightly modified version of ASQ) smearings are defined in Ref. pZ3"] . 

Straightforward application of the expansion algorithm above is theoretically 
possible but practically unfeasible: the number of monomials is huge and the 
memory requirements of the HiPPy code quickly become excessive. We get 
around this by taking advantage of the two-level structure inherent in the defi- 
nition of the action and that the intermediate reunitarisation. This allows us to 
express the partially-smeared gauge field as a member of the associated gauge 
Lie algebra, allowing us to split the derivation and subsequent application of 
the Feynman rules into two steps. 

In the first step, the Feynman rules for the outer (or "top") layer, the ASQ' 
action, are derived in the same way as before. Representing the reunitarised 
(and so far uncalculated) FAT7R smeared links by a new, Lie-algebra-valued 
gauge potential, B n : 

Ul AT7K {x) = (P U{3) o F FAT7 )[Ly = e fl M<»+5« , (37) 

we can use the HiPPy code to expand the ASQ' action in terms of as before. 
We obtain similar position-space contributions 

V r = '—^f^ Q i>(Xr;i)B^ 1 (Vr;i,l)---B^ r (v r . l . r )T r ; i 'lp(y r ; l ) (38) 

i 

that Fourier transform to give monomials of the usual form. 

To complete the derivation of the HISQ Feynman rules, we also need to 
know the expansion of B u in terms of the original gauge potential An. To 
obtain this, we write inner (or "bottom") layer, the FAT7-smeared link, as 
-?FAT7 [U ] = M — HW, where W = H and W G £7(3) (we suppress Lorentz and 
lattice site indices in the following). We again use the HiPPy code to obtain 
an expansion 

M = c[l + a n A u + a au A u A v + ..] (39) 

where c is a normalisation constant which, for simplicity in the text, we assume 
has been rescaled to unity. The notation here is, for instance, 

a^ApAv = ^2 a ^( x > v)Au(x + ^fi)A u (y + *i>) . (40) 
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Then unitarity of W implies that R = MM^ = H 2 and hence W = Rr x l 2 M 
using the expansion 

R- 1/2 = (1 + (R - l)r 1/2 = 1 - l(R ~ 1) + l(R - l) 2 + . . . (41) 

2 8 

Rearranging the result as W = exp(B), i.e. 

B = \og{W) = (W- 1) - -(W - l) 2 + ... (42) 

finally yields the desired expansion of B. These formulae are implemented in 
this form in the HiPPy code. 

Given this, we can now numerically reconstruct the HISQ Feynman rules for 



any given set of momenta from Eqn. ( 38 ) by a convolution of the ASQ' Feynman 
rules of Eqn. ( 38 ) with the expansion of B^ in terms of , summing up all the 



different ways in which the gluons going into the vertex could have come 



from the fields B^ appearing in Eqn. (38) 



In more detail, we now separately have the expansions of the ASQ' action in 
terms of the fattened gauge fields B, and of B in terms of the unfattened gauge 
field A. To obtain the correct reduced vertex, we must find all the ways that we 
can get unfattened gluons of the correct Lorentz polarisations ( "directions" ) . In 
doing this, we bear in mind that B^ contains, in principle, gauge fields A v in 
all directions and not just v = /i. Below we give an expression for the reduced 
vertex for general r. Explicit formulas for r < 3, as implemented in the HPSRC 
code, are given in Appendix |b| Using the partitions P^ as before, the reduced 
vertex for a two-level action is:: 

Yf>(p,<?;fci, mi; Mr) = 

^2 ^2 ZF,\p\(p,q;k Pl ,is 1 ;...;kp lPl ,v\p\) 

x II X £|Q|( fc Qi'^Qi;---; fc Q|Q P MQ| Q |) (43) 
QeP 

where i — 1 . . . \P\ is again the position of element Q in the ordered set P. In 
the algebra reduced vertex X, the momenta kQ i are each one of the arguments 
to the overall reduced vertex Y. As before, in the field reduced vertex Z, kp. 
refers to the summed momenta for the partition Pj . 

In Fig.[T]we represent this expansion graphically for r = 2. Solid circles rep- 
resent the Yp, r , whilst crossed and open circles represent Zp, r an d Xf,t respec- 
tively. Top-level, fattened gluons B^ are represented by dashed (brown) lines, 
whilst bottom-level, unfattened gluons A^ are shown as wavy (blue) lines. The 



two sub-diagrams represent the two partition contributions listed in Eqn. (53 1. 

As we shall discuss later, this partitioning translates naturally into blocks 
of code. For certain calculations, symmetries of the Feynman diagram will lead 
to the contributions of some of these blocks being zero. We can then improve 
the performance of the code by commenting them out in these circumstances. 
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k| i Mi k 2 , M2 k-j , p.-) k 2 , H2 




I ki + ka,^ ki Vi N^ x k2 V2 



Figure 1: A gra phic al representation of the reduced vertex Yp ir for a two- level action, with 
r = 2 as in Eqn. j53| . Solid circles represent the Yjr,., whilst crossed and open circles represent 
r and Xi? j7 . respectively. Top-level, fattened gluons are represented by dashed (brown) 
lines, whilst bottom-level, unfattened gluons Ap are shown as wavy (blue) lines. The two 
sub-diagrams represent the two partition contributions listed in Eqn. ( | 53 [ > - All momenta are 
incoming to the vertex. 



For instance, in the "tadpole" diagram in the one-loop fermion self energy, we 
can remove the term in Yp^ that is proportional to X F \ |48j . Similarly, in the 
calculation of the radiative corrections to the gauge action [33], we can suppress 
the term proportional to X F 1 3 in the "octopus" diagram. 

We can further optimise the calculation of the reduced vertices Yp, r by 
reusing terms X Fr that appear multiple times in the expressions. For instance, 
in Eqn. (54) in Appendix Ib] we can reuse X F1 (ki,[ii) and X F1 (k s , 



For tesUng it is useful to define a simpler variant of the smearings described 
on page 3 of Ref. [35] , which we call "FAT3". It is composed just of a central 

2 J c 3 — 12- 



link and adjacent 3-staples with weights c\ = i, C3 — 1 



3.2.1. Zp.r- the field reduced vertex 

The field reduced vertex essentially calculates the Feynman rule for obtaining 
r fattened gluons in particular directions fii, . . . , fj, r . The effect of fattening is 
treated separately in the X Fr . 

Zp^ r is composed of a sum of n r monomials, each representing a different 
way of obtaining the required gluons from the contours comprising the action: 



ZF,r(p, q;ki,ni;...;k r ,(i r ) = 



x + q y + ^kj ■ v nJ ] ] (44) 



where T n is a spin matrix. 

Note that if there is no fattening of the underlying gluons, this is precisely 



the reduced vertex Yp iT described in Eqn. (17 1, if we incorporate the 1/r! factor 



that was in Eqn. (14 1 



3.2.2. Xp r : the algebra reduced vertex 

If the gluons are fattened, the effect of this is contained in the reduced 
vertices Xp jr . These represent how we would choose r unfattened gluons in 
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directions fx%, . . . jj, r from a fattened link in the v direction: 




X V F (fei,//i;...;fe r ,/i r ) = - > f„exp - > k 3 -v n<j (45) 



Note that there is no spin matrix in X, so the order of the multiplication in the 
expression for Y does not matter. 

The position vectors v' n i = v n i — e„ refer to the position of the unfattened 
link relative to the start of the fattened. As such, we need to remove the half- 
link shift that was implicit in the definition of v n ^ as the centre of the link in 
the basic expansion algorithm. 

The start and end sites of the contour x and y do not appear in X , making 
X£ r=0 = lforr = 0. 



Of course, the actual values of n r , f n and {iy r } will be different in Eqns. (44 1 



and (45 1, as they come from different smearing prescriptions. 



3.3. Hardwired reunitarisation 

We described above how splitting a fermion action such as HISQ into two 
parts simplifies the generation of the Feynman rules. This allowed us to ex- 
ploit the reunitarisation of the inner FAT7R smeared links allows them to be 
expressed in terms of a new gauge potential B (suppressing the Lorentz index) . 

Nonetheless, the Feynman rules for B (as calculated in Xp^ r ) are still too 
complicated to derive for r > 3. The reason is clear if we look a tthe for- 



mula for the reunitarisation. Using Eqn. (39 1, the reunitarised field is W 



(MM') 2 M, which we can express as an element of the algebra B = log W: 



B = 


Mm + 




K = 






b^ v = 


1 / 






1 / 

2 [fi^iva 





\ ( a M l/W + al a ] + [a^ + aj^] a a ) + ^< 



(46) 



If there are n monomials in a M , there will be at least n 3 in b^ a arising from the 
term a^a^a^. There will be some compression, but this will at best only reduce 
this number by a factor of around 2. For the simpler smearing FAT3, n = 19 
and already the HiPPy code struggles to produce 6 AiyCT for FAT3R. For FAT7, 
n = 135 and direct reunitarisation to produce FAT7R links using the HiPPy 
code is out of the question. 

The alternative is to use the HiPPy code to produce the {a} and to do 



the reunitarisation in the HPSRC code by hardwiring the formulae in Eqn. (46 1. 
As before, the Xp, r algebra reduced vertices are based on expansion coefficients 
{&}. The inputs from the HiPPy code, however, implement the coefficients {a}. 
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We use these to build a set of algebra reduced vertices Wf,t and Eqn. (46 1 gives 
the relation between the reunitarised Xp r and the Wf,t' 



X F>1 (k 1 , l x 1 ) = Wp !l (k 1 , f x 1 ) 
X^ 2 (k 1 ,m;k 2 ,H2) = ^ (^,2(^1^1:^2,^2) - Wp >2 (k 2 , /i 2 ;fci, Mi)) 

^,3( fe i)Mi; fc 2,M2;fe 3 ,M3) = 

I fe^ii Mi; fe 2, ^2;k 3 ,fi 3 ) + Wp 3 (k 3 , /i 3 ; fc 2 , /i 2 ; fei,Mi)) 

-j(^F,i( fc i^i) [^ 2 (fe 2 ,M2;fe3,M3) + ^, 2 (fe 3! /i3;fe 2 ,M2)] (4?) 
+ [H^ i2 (*i,Mi;fc2,/*2) + W^ j2 (fca,A»2;*i,Mi)] ^f.iOs,^)) 

+ ^^ >1 (fel,Ml)^ 1 (fe 2 ,/i 2 )^ 1 (fe 3 ,M3) • 

4. Description of the software 

4 J. T/ie HiPPy code 

The HiPPy code is used to Taylor expand lattice actions, producing output 
files encoding the monomials making up the reduced vertices that can later be 
used by the HPSRC code to evaluate Feynman diagrams and integrals. 

The code is written in Python, for which interpreters are freely available for 
a wide range of computational platforms [50] . 

4- 1-1- Installation and compatibility 

The HiPPy code can be run on any machine which has Python version 
2.5.x installed on it. It is also expected to work with any version 2.x, but is not 
(yet) compatible with version 3.x. Only the packages supplied in a standard 
installation of Python are required. 

Installation consists of unpacking the source files in the chosen location. 

4-1-2. Testing 

A set of unit test programs are contained in subdirectory tests. They are 
run from the command line with no options and the code is expected to pass 
all of these. No input files are needed for the tests. 

Whilst the unit testing does not cover each code feature separately, those in 
test_class_f ield.py do cover most features combined. 

4-1-3. Main programs 

The HiPPy code consists of the following main programs that make use of 
the core routines described below. 

A full list of options for each of the main programs can be obtained by 
running them from the command line with option — help. 
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sample_link.py. Used to generate monomials from the algebra of the group 
for a single fattened link in the specified direction. These are used to evaluate 



X F ^ r in Eqn. (451 or W F , r in Eqn. (47) in the HPSRC code 



The fattening style is specified using — bottom_style. The output filenames 
have the form algebra_°/ s7 o i-qq* . in, where %s is replaced by the name of the 
fattening style and 7,i by the direction of the link (with value parameters. D 
mapped to zero) . Wildcard * is r copies of the character "g" (up to maximum 
parameters .R). 



For monomials that will be used in Eqn. (47 I (with hardwired reunitarisation 
carried out in the HPSRC code), the option — fld_as_alg should be specified. 
In this case 7.s is a concatenation of character "F" and the name of the fattening. 

If an undefined fattening style is specified with — bottom_style, a list of 
acceptable styles is printed. 

sample_naive.py. Used to generate monomials for the naive (or, equivalently, 
staggered) lattice Dirac action using fattened links. These are used to evaluate 



Y { F k) r in Eqn. (|35| in the HPSRC code. 



The fattening style(s) are specified using — top_style and — bottom_style, 
and should not include fattening already carried out at the algebra stage above. 
The output filenames have the form vertex_°/ s_qq* . in, where °/ s is replaced by 
the name of the top fattening style and * is r copies of the character "g" . For 
complicated fattening prescriptions, there is likely to be some filename clashes, 
so care should be exercised. 

If an undefined fattening style is specified with — top_style, a list of ac- 
ceptable styles is printed. 

sample-glue. py. Used to generate monomials for a range of simple gauge ac- 
tions, specified by a command line argument. As described above, these vertices 
have been (partially) symmetrised. 



4--1-4- Example of use 

As an example, here are the commands one would use to generate the four- 
dimensional HISQ reduced vertex files using a two-level action with hardwired 
reunitarisation and mass am — 0.2: 



python sample_link.py — bottom_style fat7 — fld_as_alg 1 
python sample_link.py — bottom_style fat7 — fld_as_alg 2 
python sample_link.py — bottom_style fat7 — fld_as_alg 3 
python sample_link.py — bottom_style fat7 — fld_as_alg 4 
python sample_naive .py — top_style asq_for_hisq 0.2 

which creates files algebra_Ff at7°/ i_qq* . in and vertex_asq_f or_hisq_qq* . in. 

The same Feynman rules can be generated using 

python sample_naive .py — top_style asq_for_hisq \ 
— bottom_style fat7r 0.2 
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but this generates far more monomials, which soon becomes prohibitive as r 
increases, as discussed above. 

For HISQ fermions, there is the option to set the mass-dependent pa- 
rameter e to zero for small am. This is done using a command line option 
— hisq_eps_zero for sample_naive .py. 

For main programs to calculate Feynman rules for other actions, e.g. NRQCD, 
interested readers should contact the authors. 



4-1-5. Core routines 

These main programs make use of the following core routines: 

parameters. py. Defines physical and numerical parameters associated with 
the expansion. If the monomials require complex-valued amplitudes, e.g. for 
mNRQCD actions, this is enabled here. 

template. py. Defines some link fattening prescriptions that can be used to 
define actions. Each fattening is a set of paths which do not need to be gauge 
covariant, so three-link Naik terms can be included, for instance. Instructions 
for adding new fattenings are included in this file. 

spin-multiplication. py. Implements the implicit representations of the spin 
algebras described in Sec. |2.6| Running this from the command line displays 



the multiplication table from Eqn. (28 1 for either Pauli or Dirac spin matrices. 



wilson. py. The main expansion engine that converts a collection of Wilson 
lines into a Taylor expansion in the form of reduced vertices built from mono- 
mials. The collection of Wilson lines is defined by 

thepath = [ [cl , c2 , \dots] , [pi ,p2 , \dots] ] 

where cl, c2 etc. are the amplitude of each contour, the path of which is defined 
by a set of signed integers. For instance, a plaquette would be pl= [1,2,-1, -2] . 

As an example of a complete path, a two-dimensional naive fermion action 
would be defined by a path: 

thepath = [ [0. 5, -0.5, 0.5,-0. 5, m] , [ [1] , [-1] , [2] , [-2] , [] ] ] 

The expansion is not symbolic, so the mass m must take a specific numerical 
value. 

It is assumed that all links are of the same, fattened style. The fattening 
can contain any number of iterated styles defined in template .py, given as a 
list top_style interpreted from right to left. 

At the lowest level of this iteration, a further single level of fattening can 
be specified using bottom_style. Uniquely, a reunitarisation by projection onto 
U (N) can be imposed after this fattening by appending character "r" to the end 
of the name of the style (again, chosen from the definitions in template .py). 
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cl ass .entity .py. Each monomial in the reduced vertex is encoded as an in- 
stance of this class, a description of which can be found in Sec. 4 of Ref. [T7] . 
Translation invariance is exploited so that all paths are assumed to start at the 
origin, x = 0. If this is not the case (e.g. for open boundary conditions in 
the temporal direction for Schrodinger functional calculations [3EJIMI), then an 
attribute start_site would need to be introduced to this class and appropri- 
ate changes made to ensure gauge covariance in the multiplication and in other 
places. 

class_field.py. The collection of monomials defining the reduced vertices re- 
turned by wilson. path.0 are stored as a dictionary in an instance of this class. 
A fuller description can be found in Sec. 4 of Ref. [i~7] . 

class_algebra.py. Very similar in form to class_f ield.py, this class contains 



the reduced vertices associated with the algebra in Eqn. (42 1 



mathfns .py. contains a collection of methods not tied to any particular class. 

4.2. The HPSRC code 

The HPSRC code uses the input files from the HiPPy code to create Feyn- 
man rules, and combines these to evaluate Feynman integrals by exact mode 
summation or statistical estimation using the Vegas algorithm [14j [15] . 

4-2.1. Installation and compatibility 

The HPSRC code requires a standards-compliant Fortran95 compiler with 
support for minimal CPP directives (chiefly #ifdef and #ifndef). The com- 
mand make is optional but useful. An appropriate MPI wrapper for the compiler 
will be needed for creating a parallel version of the executable. No additional 
libraries are required. 

The code is known to compile and execute correctly with the following 
compilers: Intel ifort, Portland pgf90, NAG f95 and GNU gfortran on Unix- 
based systems (the first two also with MPI on parallel machines) and Silver- 
frost/Salford ftn95 and Lahey-Fujitsu lf95 on Microsoft Windows systems using 
cygwin. 

To install the code, the source should be unpacked in a suitable directory 
and the file Makef ile_details (which contains machine-dependent information) 
edited to point to the correct compiler. 

It is a good idea to execute make clean before compiling the code. A given 
target this_executable can be compiled using the command 

make this_executable 

Use the command make to see a list of possible executables, or examine Makefile. 



23 



4-2.2. HPSRC filename conventions 

Main programs are named main_* . F90 and compile to form executables with 
the same filestem and extension .$(EXE) as specified by Makef ile_details. 
Typically $(EXE) = x, but certain compilers expect $(EXE) = exe. 

A Fortran module this is typically contained in a file named mocLthis . F90. 

Input files containing information that should be specified at compile time 
are named with extension .i. The command make clean must be executed 
before recompiling if any of these files are changed. These files typically con- 
tain unavoidable CPP directives and the number of these has been kept to a 
minimum. 

Input files that are read at runtime arc have names with extension . in. Most 
input files have a filename that reflects that of the module that reads them. 

4.2.3. Testing 

The code has been rigorously tested as a whole in a number of scientific cal- 
culations, being compared with identical calculations carried out using different 
programs. 

In addition, various manual and automatic tests are routinely carried out to 
verify the continued correctness of the code. The manual tests are located in 
subdirectory tests/. 

Automatic differentiation. Without a tcmplating option in Fortran, each vertex 
is effectively coded twice: once as a complex-valued object and once as a taylor- 
valued object that includes the automatic derivatives. 

The code exploits this by performing a finite difference of the first to compare 
with the derivatives included in the second. To avoid multiplication of errors in 
higher-order difference operators, the tests are carried out as follows. First, the 
non-taylor version is compared with the zero-order derivative from the taylor 
version. If this agrees, a first-order, symmetric difference is made of the zero- 
order derivative in the taylor in direction r. We check that, for all r, this 
agrees (within a specified tolerance) with the analytic first order derivative in 
the taylor object. If so, we form first-order differences of the analytic first-order 
derivatives and compare these with the analytic second order derivatives and so 
on. 

These comparisons are done at runtime when the vertex modules are first 
initialised. 

This test checks not only that the automatic differentiation (and associated 
overloading of operators) is working correctly, but also that there are no coding 
inconsistencies between the taylor and non-taylor versions of the vertices. It 
cannot, of course, detect algorithmic errors that have been consistently coded 
in both versions. An additional, lower level test of the TaylUR package is 
provided by test_taylors .x. 

test_compare_vertices .x. Can be used to check that the Feynman rules are 
the same for two different implementations of the same quark action. This pro- 
vides a rigorous check of the algorithms in mod_vertex_qq* .F90. Having com- 
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piled test_compare_vertices . x, a Python script test_compare_vertices . py 

can be executed to complete various independent checks of the two-level action 
decomposition, the hardwired reunitarisation and the summand/f actor division 
of the action. The last is based on an NRQCD action and also provides a strong 
test of the implementation of the Pauli spin algebra. 

Tests of Vegas. The Vegas code is automatically tested at runtime, evaluat- 
ing the area under a two-dimensional, normalised Gaussian as well as a contour 
integral. These tests can also be run using test_vegas . x. 



Lower level HPSRC tests. test_spinors . x tests type spinor defined in Sec. 2.6 
In addition, test_pauli_dirac_spinors . x tests the embedding of Pauli two- 
spinors into Dirac 4-spinors. test_print_vertices .x prints the values of the 
fermionic vertices for a random set of momenta. The same random number seed 
is used each time, so running the code for different vertex_qq_composite . in 
on the same machine will yield outputs that can be directly compared. 

4-2.4- Main programs 

To illustrate the use of the HPSRC code, we provide three sample programs 
with Makefile targets: 

quark_sigma. x. Calculates the one- loop renormalisation of ASQTAD improved 
staggered fermions. To do this, it evaluates the self-energy Feynman diagrams 
and their derivatives. 

nrqcd_sigma.x. Performs a similar calculation for heavy NRQCD fermions. 

nflwimp . x. Calculates the HISQ fermion contribution to the one- loop radiative 
corrections to the Luscher-Weisz gauge action, as per Refs. [M| [35] . 

4-2.5. Core routines 

This subsection describes modules that are needed to use the vertex modules. 
If MPI is to be used, the file use_mpi . i must be changed prior to compilation. 
Similarly, if monomials require complex amplitudes, e.g. for mNRQCD calcula- 
tions, this is enabled in use_complex_amplitudes . i. 

mod_num_parm.F90. Numerical parameters that are used by most of the other 
modules. Most parameters in this module will not need to be changed. It also 
optionally reads from file paths . in the pathnames of the directories holding 
the vertex and algebra files generated by the HiPPy code. 
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mod_phys_parm. F90. Physical parameters used in the calculation. The code 
has only been tested for number of dimensions Ndir = 4 and colours Ncol = 3. 
Other choices will require careful testing before use. Values for many parameters 
can be changed at runtime by editing file phys_parm. in Twisted boundary con- 
ditions are selected by setting the number of twisted directions ntwistecLdirs 
to 2, 3 or 4. Value gives periodic boundary conditions. For finite lattices, 
the lattice size is specified in latt_size(0:Ndir-l), with the first component 
the temporal one. If momenta are to be squashed — > fc M — a (tt sin(fe /J ) [18], 
the squash factors in each direction are specified in mom_squash(0:Ndir-l). 
The anisotropy metric anis is only supported currently for gluonic vertices. 

mod_momenta.F90. Defines type mom which encodes momenta. Each instance 
contains the momentum components, the twist vector (only relevant for twisted 
boundary conditions) and a route variable that is used in automatic differenti- 
ation. 



mod_taylors .F90. Defines type taylor, as described in Sec. 2.5 



modspinors . F90. Defines types spinor and tayl_spinor, as described in 
Sec.EU 



mod_matrices . F90. Defines types mat4x4, used to explicitly represent objects 
such as the gluon propagator in Lorentz index space. A corresponding taylor- 
valued type tayl_mat4x4 is also defined. 

mod_colour .F90. Calculates the Clebsch-Gordan factors for both gluonic and 
fermionic vertices, for periodic and twisted boundary conditions. Colour factors 
for periodic boundary condition are pre-computed and stored in the compile- 
time files tr_cols_N. i and mat_cols_N. i. For twisted boundary conditions, the 
Clebsch-Gordan factors are computed as needed. 

mod_mod_mpi .F90. The code can either be run as a scalar code on a single 
processor or as a parallel code using MPI. To achieve this with the minimum 
of compile-time code modifications (either by hand or using CPP directives), 
mod_modjmpi . F90 contains a set of dummy MPI subroutines. Whether these 
dummy routines are used or not is controlled by the single CPP directive in file 
use_mpi.i. Note that mod_mod_mpi .F90 is the only file that reads use_mpi.i. 
Only the MPI commands we use are included in mod_mod_mpi . F90. We do not 
know of any externally maintained library of dummy MPI commands for use in 
such situations. If this were to exist, it could easily be used to replace this file. 

Throughout the code we are (or try to be) careful that input/output files 
are only opened by one processor, with information then shared using MPI 
commands. 
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mod_vegasrun . F90. Provides an implementation of the Vegas algorithm, based 
on the original code of Peter Lepage [HI [T5] , who also helped convert the code 
to a parallel version using MPI. On parallel machines, the processors are divided 
into farms, each of which does an independent Vegas estimate of the integral. 
Within the farm, all processors carry out the same calculation, interacting only 
to share the work of evaluating the integrand at each of the Monte Carlo points. 
The number of processors per farm is controlled in file vegas_parm. in, and 
must be a factor of the number of processors. Farming has the advantage of 
avoiding potentially slow global communications on large clusters. 

An improved parallel implementation is currently being prepared as part of 
a DEISA DECI-funded project [51, . 

4-2.6. Gluonic vertex modules 

Feynman rules Vc,r for the gluonic vertices are implemented in 
mod_vertex_* . F90, where * is r (the number of gluons) copies of character "g" . 
Currently vertices for r < 5 have been implemented. 

The modules for r > 1 are called vertex_* and have a very similar structure. 
The top level routine is vert_*(k, lorentz, colour), which takes as arguments 
the momenta, lorentz and colour indices and returns the complete symmetrised 
Feynman rule for that vertex as a complex number. The colour array is ig- 
nored if we are using twisted boundary conditions. The module assumes that 
the monomials written by the HiPPy code have been partially symmetrised as 
described above. The remaining symmetrisation over the distinct cosets is car- 
ried out automatically in the top level routine. The reduced vertex is calculated 
from the monomials in yvertex_*(). 

This code structure is repeated using instances of type taylor instead of 
complex numbers, providing analytic derivatives of the Feynman rules. The top 
level routine in this case is called taylor_vert_* () . As detailed above, coding 
the algorithm twice in this way is exploited in the runtime testing of the code. 

The remainder of the mod_vertex_qq* . F90 files is taken up with initialisation 
code (which reads in all the vertex files at runtime) and testing routines. 

Only mod_vertex_gg.F90 differs from this overall structure. This calculates 
the gluon propagator. In this case, the top-level routine is 
gluon_prop(k, colour) . The propagator is simultaneously calculated for all 
Lorentz index combinations, packaged as a D x D matrix implemented as derived 
type mat4x4. 

The choice of gauge is controlled by variable Galpha, specified in runtime 
input file vertex_gg_parm. in. Value 1 corresponds to Feynman gauge and 
to Landau gauge. As detailed in Ref. [28 , the inverse propagator does not exist 
in Landau gauge, so an intermediate gauge parameter Gbeta should be used. 
gauge_f amily = 'landau' is the usual choice, where the gauge correction is 
based on the four-vector k. Changing this to 'coulomb' uses just the spacelike 
components instead, but this option is not fully tested. 

As evaluation of the gluon propagator can be an expensive part of a per- 
turbative calculation, the propagators for various gluon actions are hardwired 
in mod_vertex_gg.F90. Their use is specified by changing variable action 
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in vertex_gg_parm. in from 'python' to appropriate other character strings 
which are listed in subroutine vertex_gg_params () . 

A gluon mass can be added to the propagator. Its square is gluon_mass2, 
specified in vertex_gg_parm. in. 

Again, this code is repeated using automatic derivatives packaged instead 
in a taylor-valued matrix of type tayl_mat4x4. The remainder of the module 
contains runtime tests and initialisation routines. 

Associated with the gluon modules are modules calculating the Feynman 
rules associated with the Fadeev-Popov ghost fields and the Haar measure. 
These arc implemented as hardwired formulae in mod_vertex_hhstar .F90 (for 
r < 2) and in mod_vertex_meas_gg.F90 (only for r = 2). 



4-2.7. Fermion vertex modules 

Feynman rules Vp> for the fermionic vertices (the unsymmetrised version 
of Eqn. (15l) are calculated in the files mod_vertex_qq* . F90, where * is r (the 
number of gluons) copies of character "g" . Currently vertices for r < 3 have 
been implemented. The modules for r > 1 are called vertex.qq* and have a 
very similar structure. The top level routine is vert_qq*(k,lorentz, colour), 
which takes as arguments the momenta, lorentz and colour indices and returns 
the complete Feynman rule for that vertex as an instance of type spinor. The 
convention is that k(l:r+2) is (/q,p, fei, ■ ■ ■ , k r /). The colour array has the 
same ordering, but is ignored if we are using twisted boundary conditions. 

Internally, this function calculates the reduced vertex using yvertex_qq* () 
and multiplies on the appropriate Clebsch-Gordon factor, calculated using 
mod_colour . F90. 



yvertex_qq* () . Implements the summands and factors decomposition described 
in Sec. 13.11 The reduced vertices for each factor are calculated in 
yvertex_qq*_partial(). The sum over partitions X)p<ep<>-) naturally translates 
into a block structure for this routine, the order of which follows that of the 
expressions in Appendix |A"| 

yv er t ex -qq* -partial () . Implements the two- level field/algebra algorithm as 



per Section 3.2 The "upper-level" , fattened, field vertices Zp r are calculated in 
yvertex_qq*_partialjnoalg() , combining monomials that have been read into 
the module at runtime. The "lower-level", algebra vertices Xp, r are calculated 
in yvertex_qq*_algebra() . Once more, the sum over partitions translates into 
a block structure for the code, and the order reflects the expressions in Ap- 
pendix [B] As discussed earlier, for certain Feynman diagrams we can comment 
out some blocks on symmetry grounds. 

yvertex_qq*-algebra() . Calculates Xp, r in terms of the monomial-calculated 

There is the option to 



Wp r from y_qq*_alg_basic(), as detailed in Sec. 3.3 
apply no reunitarisation (setting Xp. r = Wf,t)i or t 
described in the text. Additional projection methods, such as "stouting 
could also be implemented here. 



re hardwired projection 
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This code structure is repeated using instances of type tayl_spinor instead 
of type spinor, providing analytic derivatives of the Feynman rules. The top 
level routine in this case is called taylor_vert_qqg* () . As detailed above, 
coding the algorithm twice in this way is exploited in the runtime testing of the 
code. 

The remainder of the mod_vertex_qq* . F90 files is taken up with initialisation 
code (which reads in all the vertex files at runtime) and testing routines. 

Only mod_vertex_qq. F90 differs from this overall structure. This calculates 
the fermion propagator. In this case, the top-level routine is 
quark_prop(k, colour) where the momentum specified is q (which flows in the 
direction of the fermion arrow). inv_quark_prop ( ) calculates the two-point 
function for the fcrmions and implements the summands and factors decompo- 
sition described in Sec. 13.11 The reduced vertex for each factor is calculated 
in inv_quark_prop_partial () . This is the lowest-level routine, as there is no 
algebra contribution for r = 0. These routines are followed in the code by type 
tayl_spinor versions and by initialisation and tests. 

4-2.8. Runtime specification of Feynman rules 

The HPSRC code allows for the use of multiple fermion types, each with 
their own Feynman rules constructed in their own way. The details of all of this 
are specified at runtime by the file vertex_qq_composite . in, which is read by 
all the mod_vertex_qq* . F90 files. The input takes the form of a NAMELIST. 
A typical example is given in Appendix |C| 

5. Conclusions 

The accuracy of simulations of lattice QCD (and other field theories) has 
been markedly improved by the use of Symanzik- and radiatively-improved ac- 
tions and measurement operators. In addition to the calculation of the radiative 
corrections, this improvement programme also requires a concomitant effort in 
calculating associated renormalisation parameters that are vital to the contin- 
uum and chiral extrapolations. Lattice perturbation theory is an essential tool 
in these calculations. For some quantitities we can use other techniques, in 
particular non-perturbative renormalisation (NPR) [TJ high-/3 simulations 
[TU1 [T2] and stochastic techniques [13 . In cases where these can be used, lattice 
perturbation theory provides a valuable check of their results and, in the case of 
high-/3 simulations, important constraints on the coefficients of the low powers 
of a s . When there is complicated mixing of operators, which happens in a lot 
of physically relevant measurements, statistical errors lead to the failure of the 
NPR and high-/? techniques and lattice perturbation theory is the only alter- 
native. Stochastic perturbation theory is also very expensive for theories with 
dynamical fermions. 

Lattice perturbation theory for improved actions is complicated: not only 
does it include many vertices not present in the continuum, the Feynman rules 
also contain a very large number of terms. These complications vanish, of 
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course, in the continuum limit as the lattice spacing a tends to zero. The point 
is, however, that lattice perturbation theory is used to correct for the missing 
momentum modes on a discrete lattice, so the presence of these complications 
are central to its utility. Many of the complicating terms violate Lorentz sym- 
metry and are mathematically complicated. This makes analytic evaluation of 
Feynman integrals impossible in almost all cases. It is therefore crucial that 
we develop a robust computational method for deriving and then using these 
Feynman rules in perturbation theory. This method must be efficient as well as 
accurate. 

In this paper we have described such a method. We have used a two-pronged 
approach, with separate programs to expand the action and derive the Feynman 
rules, and then to use these rules to calculate Feynman integrals. In both 
cases, we have developed efficient algorithms to minimise the computational 
expense. Indeed, without these efficient algorithms we find the calculations 
to be impossible for many realistic action choices. We have also presented 
a full, working implementation of the algorithm in the form of the HiPPy 
and HPSRC codes. The programming languages used for the implementations 
are Python and Fortran95, respectively. Python offers good list-handling and 
object-orientation features that are suited to the expansion of the action to 
derive the Feynman rules. Fortran95 is numerically efficient, an over- riding 
concern for the evaluation (or estimation) of expensive Feynman integrals. 

We stress that our method is gcncrically applicable. In the text, we have 
focused on the features of our method that make possible calculations using 
realistic NRQCD and HISQ fermion actions. This bias merely reflects the prob- 
lems for which we have used the method most so far. We have considered other 
actions, both by expanding them using the HiPPy code and by hardwiring the 
hand-derived Feynman rules in the HPSRC vertex modules. 

Particular strengths of our approach are the ability to deal with complicated 
colour and spin structures in actions and the breaking down of actions into sim- 
pler sub-structures for faster evaluation. We have also provided a particularly 
simple way of describing actions that reduces the scope for transcription errors. 
The comprehensive suite of tests in the code also give confidence in the accuracy 
of the basic components of both the HiPPy and HPSRC code suites. 

As well as extending the basic functionality of the codes, e.g. to describe 
actions with multiply nested levels of link fattening, the codes are easily ex- 
tended to a wide range of related problems in lattice simulation, as has already 
been done for lattice chiral perturbation theory |37j and Schrodinger functional 
calculations [35J ■ Another area in which the code might be used is in cal- 
culations where the lattice gauge field is split into a fluctuating, quantum piece 
and an external background field |53[ 154] . The authors are happy to provide 
further details on the implementation of this, as well as additional advice on 
the use of the HiPPy and HPsrc codes. 
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n 1 n 2 n 3 



Figure 2: A graphical representation of the reduced vertex Yp r for r = 3 for a single summand 
of an action that is composed of a number of factors, as described mathematically in Eqn. |51| 
and implemented in the HPSRC code. The four terms in Eqn. ( |51| are separately shown, with 
"1,2,3 showing the number of the factors from which each of the gluons are drawn. Open 
circles denote vertices associated with individual factors. Factors contributing no gluons are 
show n as "— • ■ ■ — ". Momentum flow in the vertices is not shown, but can be deduced from 
Eqn. | |5l} . 
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7. Appendices 

A. Explicit expressions for factors 

Using Table [l] we here give explicit expressions for Eqn. d35| for Y Fr for a 
summand in the action with N factors, for r = ... 3 (the range of r imple- 
mented in the HPSRC code): 

JV 

Y F , (p,q) = l[Y$(p,q) (48) 
fc=i 
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N /m-l \ 
Y FA (p,q;k 1 , f i 1 )= J2 II y^(P,-p))y^i\p,Q;k 1 ,fi 1 ) 
m=l \ fc=l / 

X ( ft ^S(-9.«)J ( 49 ) 

\fc=7li + l / 

yF,2(p,g;fci,Mi;*!2,/i2) = X! II ( p > _p ) J 

i<m<JV \ fe=i / 
x Y^ 2 1 \p,q;k 1 ,fi 1 ;k 2 ,^2) f II y i!o( _ 9>9) I 

\fe=ni+l / 

e ( n ^(f'-PH^iP'-p-^^^o 

l<rii<n 2 <A r V fe=l / 

( n Y^ip+ku-p-kM 

\k=m+l / 

xy^ 1 2) (p + fe 1 ,g;fe 2 , Al2 ) ( [] y i3(-9^)) (50) 

\fe=n 2 +l / 
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E J^ Y F^P>-P)) Y F^HP'9-,ki,fJ,i-,k2,lJ.2-,k s ,iJ, 3 ) 



i<m<N \ k=l / 

/ N > 

x n y $(-q,<i) 

frt\ — X 



l<m<n 2 <N \ k=l / 



x ( II ^o(-g- fc 3,g + fc3)) 

\fc=m+i / 
<Y^\-q-k a ,q;k 3 ,n 3 ) ( J] Y%(-q,q)\ 

\k=n 2 +l / 

E ( fi Y F k o(p, -p) J ^ V -p - 



l<ni<Ti 2 <A r \ fc=l 

/ "2-1 

x 



n ^i3(p+ fe i'-p- fe i) 

\fc=rii + l / 

F^ 2 2) (p+fe l! q;fe 2 ,M2;fe3,/i3) ( [] F i3(-<7,<?)) 

\fc=n 2 + l / 

e ( n ^(p.-pH^fe-p-fe^i.ft) 

:n2<"3<A r V fc=l / 



l<rai<n 2 <Ti3<Af \ fc= 
/ "2-1 

X 



n y^ip+k^-p-k!) 



ri3 — 1 

X$'(P+*1,? + *3I*J,/U) ( II Y^(-q-k 3 ,q + k 3 )\ 

\fe=n 2 +l / 

xY^\-q-k 3 ,q;k 3 ^ 3 )( ]J Y$(-<?,<?) ) (5 

\fe=n 3 +l / 

The reduced vertex from the n th is denoted YpJ . The formula for r = 3 
illustrated graphically in Fig. [2] 

B. Explicit expressions for two-level actions 

To make things more concrete, we give explicit formulae for the r = 1, r = 



and r = 3 reduced vertices for a two- level action as defined in Eqn. (43 1. This 



33 



the range of r implemented in the HPSRC code. Using the partitions in Table [T] 
we obtain: 

Y Ft x(p,q;ki,m) = z f,i(p, q; fei, vijXp^kx, m) (52) 

YfMp> 95 fe i> Mi; fe 2, M2) = J! ^-F.iCPi 9; fc i + fc 2, u^X^ki, m; k 2 ,n 2 ) 

+ 51 ^,2(p 5 9;feij iy i;fe2,^2)^i(fei,Mi)^i( fe 2,M2) (53) 

y f ,3(p, q; k 1 , fii; k 2 , k 3 , /.i 3 ) = 

^2 Z FA (p, q; fei + k 2 + k 3 , Vi )X F 1 3 (fei , // x ; fc 2 , M2; ^3, M3) 
^1 

+ ^ Zp : 2{p,q;k 1 + k2,vi;k 3 ,iy 2 )X F 1 2 (k 1 ,n 1 ;k2,fi2)X F y(k 3 ,fi 3 ) 
+ Z F: 2{p,q\k 1 ^ 1 ;k 2 + k 3 ^2)X F 1 1 (k 1 ,^ 1 )X F 2 2 (k 2 ,^2;k 3 ,fi 3 ) 

+ z ^s(p, q; fei, fi; fe2, fa; fea, ^3)^i(fei, Mi)^j?i( fc 2, ^2)^1(^3, Ma) 

(54) 

C. Example of vertex_qq_composite . in 

Here we provide an explicit example of the runtime file used to specify the 
Feynman rules used by the HPSRC code. 

The calculation will use two types of quarks: relativistic HISQ and heavy 
NRQCD. It is assumed that the appropriate vertex and, where relevant, algebra 
files have been precomputed using the HiPPy code and stored in the locations 
specified in the HPSRC file paths . in. 

&vertex_qq_composite 

no_quark_types = 2 
! First quark type: HISQ 

! All RH indices set to 1 to denote first quark type 

summands(l) = 1 

factors (1:1,1) = 1 

summand_amps ( 1 : 1 , 1 ) = l.dO 

opnamed : 1 , 1 , 1) = "asq_f or_hisq_" 

algname(l) = "Ffat7" 
reunit_to_apply_to_alg(l) = "project" 
! Second quark type: NRQCD 
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! All RH indices set to 2 to denote second quark type 
summands (2) = 2 
factors (1:2, 2) = 0,4 
summand_amps(l:2,2) = I.d0,-l.d0 

opname(l:4,2,2) = "nrqcd_A" , "nrqcd_B" , "nrqcd_C" , "nrqcd_A" 
algname(2) = "simple" 
reunit_to_apply_to_alg(2) = "none" 

! N.B. no opname ( : , 1 , 2) definition because first summand has no 

! HiPPy inputs because it is a constant 

/ 

That there are two quark types is specified in the NAMELIST using 
no_quark_types. Within the HPSRC code we control which set of Feynman 
rules we use by setting variable quark_type to 1 or 2 prior to calling 
quark_prop() or vert_qqg*(). 

The first fermions are relativistic HISQ quarks. There is no splitting of the 
action into summands and factors, so summands (qt)=l and 
factors (1 : summands (qt) ,qt)=l. Similarly, the only summand has amplitude 
summand_amps (1 : summands (qt) ,qt)=l . OdO. The files used to construct Zp r 
for each factor for each summand are assumed to have been generated by 
the HiPPy code and stored in files named vertex_°/ sqq* . in where %s is re- 
placed by appropriate text for each factor for each summand as specified in 
opname(ft,sm,qt), where ft runs from 1 to f actors(sm,qt) and sm runs 
from 1 to summands (qt) . In this case, the only vertex files are those with °/ s 
replaced by asq_f or hisq_. 

The name of the algebra file is given by algname (qt) (note there is no trailing 
underscore in this case). We specify that hardwired projection is required to re- 



late Wf.t in Eqn. (47 1 to Xp tT in Eqn. (45 I using reunit_to_apply_to_alg(qt) 



The HPSRC code currently assumes that the same algebra files are used for all 
factors and summands of a particular quark type. Again, it is assumed that the 
files defining these have been pre-produced using the HiPPy code and stored in 
files named algebra_y o s7 i-qq* . in. °/ s is replaced by algname (qt) and 7,i runs 
from to D — 1, where D is the number of dimensions. 

In this explanation, we have used qt to denote the quark type, and other 
intuitive labels to specify array sizes. Unfortunately, NAMELISTs only support 
numbered array entries. 

The second block of the NAMELIST defines an NRQCD action of the heuris- 
tic form ^(1 — ABCAjip. There are two summands, the first with factors 
(which gives default answer 1) and the second with 4. The minus sign in front 
of the second summand is specified in summand_amps. The first summand has 
no factors, so opname is not specified. The second summand has 4 factors, and 
the filenames are specified. Note that the first and last filenames are the same. 
The algebra is specified as before, and no hardwired reunitarisation is required. 

In both cases, the correct spin algebra is specified in the headers of the input 
files written by the HiPPy code. 
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